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Introduction 

Astronomy is an area of applied physics in which unusually 
beautiful objects challenge our imaginations to explain observed 
phenomena in terms of known laws of physics. It is a field that 
has stimulated the development of lots of physical laws and lots 
of mathematical and computational methods. Many statistical 
and numerical methods may be traced to 19th century efforts to 
determine orbits of solar system objects such as comets and asteroids. 
Numerical calculations were to be done by hand, with pencil and 
paper. Methods were needed that would guarantee a reliable orbit 
with minimal effort, based on imperfect observations that were 
not evenly spaced in time. In another area, some well-known and 
powerful stability techniques were developed to study the classical 3- 
body problem by mathematicians like Poincare. The stability of the 
solar system is a problem that continues to challenge mathematicians. 
The story of Newton, apples, and the moon’s orbit is well known to 
all of you. 

Outstanding scientists were drawn to astronomical problems by 
their charm and beauty, not to deny practical usefulness, in an earlier 
era. The opening paragraphs of Maxwell’s Adams Prize Essay (1859) 
on the structure and stability of Saturn’s rings eloquently appeal to 
the charm of trying to understand a natural phenomenon for which 
no practical use could be imagined. He did that in an era when 
astronomy was the quintessence of “practical'’ science for navigation 
at sea. Techniques worked out to attack astronomical problems are 
useful in many other contexts, and they promise to remain useful for 
years to come. 

Charm and beauty characterize many problems of present-day 
astronomy as well. As in earlier eras, astronomical problems drive 
us to consider phenomena in extraordinary physical conditions. 
Physical conditions in astronomical objects take on far more extreme 
values than can be mimicked in any terrestrial laboratory. Densities, 
for example, range from less than one atom per cubic centimeter 
(p ~ 2 x 10 -24 gm/cm 3 ) in interstellar space to nuclear densities 
(p 10 15 — 10 16 gm/cm 3 ) in neutron stars, rising a bit higher if 
the star is on the verge of collapsing to a black hole. Temperatures 
range from around 100 ° K in molecular clouds to 10 10 °K at the 
centers of pre-supernova stars. Magnetic fields can be found that 
range from a microgauss in interstellar space to 10 12 gauss at the 
surface of a neutron star. Magnetic fields are important because the 
interstellar medium tends to be a highly ionized plasma, thus to have 
high conductivity. 

Boundary conditions are much cleaner in astronomical problems 
than in laboratory experiments as well. Astrophysical plasmas 
(e.g., the solar wind, solar atmosphere, interstellar medium) are not 
contaminated by boundaries as are laboratory plasmas. Problems 
in which laboratory boundary conditions differ from real-world 
situations are familiar to people working on aircraft design as well. 

Astronomers have been in supercomputing for the past century 
and a half as well, if we define supercomputing to mean numerical 
computation on an unprecedented scale for its time. Imagine the 
numerical work in computing ephemerides for navigational purposes 
in the 19th and first half of the 20th centuries. An astronomical work 
still stands as a landmark from the days when it was first done with 
punched-card calculators: the “Coordinates of the 5 Outer Planets 
1653-2060,” by Eckert, Brouwer, and Clemence (1951). A sign of the 
times: five years ago, amateur astronomers took it as a challenge to 
repeat that calculation on PC’s like the TRS-80. Fortunately, they 
could refer to the original work for methods and for numerical values 
that could be used to check their calculations. 

Three features drive present research efforts that impel as- 
tronomers to seek the computational power of supercomputers. 


(1) A complex interplay of many physical processes. Interactions are 
almost always nonlinear. General physical situations must be taken 
into account to fit observations. 

(2) Higher quality observational material discloses details of flow and 
emission patterns that require complex models. Several factors con- 
tribute. (a) Improved (multi- element, low-noise) detectors, (b) Ex- 
tension to new wavelength bands (radio, infrared, ultrarviolet, x-ray, 
gamma rays), with (c) high resolution imaging and (d) spectroscopic 
resolution as well. These features come from excellent new instru- 
ments such as the VLA and VLBI in radio astronomy, from KAO and 
IRAS and the future SOFIA and SIRTF in infrared astronomy, from 
IUE and the Hubble Space Telescope in visible and ultra-violet wave- 
lengths, and a variety of satellites for x- and gamma rays. European 
and Japanese satellites Eire coming on line as well, and we saw coop- 
erative programs with them and with the Soviets in the recent Halley 
missions. Scientific use of ail these superb instruments requires source 
modeling in at least as great detail as the instruments can yield. 

(3) Reduction of raw observational data to a form comprehensible 
by people. Instruments like the VLA (a 27-element interferometer or 
synthetic aperture) produce what is essentially a Fourier transform 
of the source brightness distribution (in two spatial dimensions for 
each spectral passband for possibly as many as 256 spectral bands), 
corrupted by sky noise, by receiver and antenna noise, by a “dirty” 
raw synthesized beam, and messed up by the atmosphere. Other in- 
struments produce data encoded differently, but the data reduction 
problem is much the same. High data rates exacerbate the problem. 
These newer instruments produce a lot of data in a short time, then 
they switch to a new source and produce a lot more data. 

We see again the pattern that techniques developed for astro- 
nomical use have uses in other areas: Modern oil field exploration 
uses techniques roughly equivalent to turning the VLA upside down 
so its antennas detect signals coming from inside the earth. Some of 
the astronomical imaging methods such as CLEAN, self-calibration, 
hybrid mapping, etc., are finding use there. Maximum entropy 
restoration methods are being developed actively in both camps, as 
they are elsewhere. 

Image restoration in radio astronomy is truly impressive. Within 
the past couple of years, a map of the radio source Cygnus-A has 
been produced with 1000 or more pixels along an edge, and with a 
dynamic range (ratio of brightest spot to faintest detail recorded) of 
about 10,000 to 1. This is in the presence of sidelobes that can range 
to 10% for a “dirty beam.” It is difficult to look at that map without 
a tingle of admiration and excitement. 


Current Computations 

Computation has had an impact on nearly every area of 
astrophysics. Many of the problems of astrophysical interest 
have their origin in trying to understand the structure and aging 
(evolution) of stars. Particularly interesting computational problems 
arise with studies of the final death- throes of a massive star -a 
supernova outburst. A core like a neutron star (density around 
10 15 gm/cm 3 , temperature about 100 Million degrees K, size about 
10 miles across, about 1 solar mass or 10 33 gm) is thought to pass 
a stability limit in its equation of state that leads to a collapse. A 
shock is sent outward into the surrounding envelope. Time scales are 
milliseconds (sound speed is nearly the velocity of light). Nuclear 
burning in the shock front transmutes most nuclei to the doubly 
magic /Vt"’ 6 (28 protons and 28 neutrons). So many neutrinos are 
produced, and they carry such a large fraction of the energy, that 
neutrino opacity limits energy transport. Things become particularly 
interesting when the shock encounters the innermost nuclear-burning 
shell of the pre-supernova star where silicon is being burnt. The 
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upshot is that a burst of energy (10 54 ergs) is produced that makes 
a star as bright as an entire galaxy (10 11 solar luminosity or about 
10 44 ergs/sec) that takes weeks or months to expand and cool, as the 
outburst gradually fades. 

This is but one of the fascinating problems astrophysicists are 
tackling computationally today. Star formation, stellar evolution, 
and physics of the interstellar medium are others. Common threads 
run through the physics of all these situations. That physics 
includes (1) Radiative transfer (2) Transient fluid flows, and (3) 
Radiation hydrodynamics. This last term refers to situations in 
which radiation pressure dominates mass motions. Actual physical 
situations involve some cases where radiation dominates, others 
where fluid flows dominate, and transition cases. The fluid flows 
influence the radiation field and radiative transport as well. Strong 
shocks frequently occur in examples of astrophysical interest. The 
state of the art can cope with these cases, and nuclear burning is 
included in modern codes. Nuclear burning provides an energy 
source (and a sink in some cases). Nuclear reaction networks bring 
in several hundred nuclear species with vastly different lifetimes, so 
the problems are very “stiff.” Magnetic fields have so far not been 
successfully included, so we cannot yet handle hydromagnetic cases 
with radiation and nuclear burning. But that is coming. It should be 
here in a couple of years. 

Numerical relativity is another area receiving a lot of attention. 
Relativity here means general relativity -well known for numerical 
intractability because the non linearities affect the geometry of 
spacetime. Black holes were made to collide a few years ago (lots 
of gravitational radiation emitted). Star clusters going relativistic 
are just now being studied-but even that with restrictive symmetries 
because of numerical and conceptual difficulties. A star cluster goes 
relativistic when its central density and size become great enough to 
be near the Schwarzschild radius for a black hole. 

We see once more that new methods have to be developed to 
handle extreme physical conditions that first come to our attention in 
astrophysical situations. Often these methods can help with problems 
in other scientific and technical areas, even if conditions are not so 
extreme. Algorithms must be robust and dependable. 


Dynamics of Galaxies 

Let me turn to another particularly beautiful area of computa- 
tional sis trophy sics, the area of my own research, the dynamics of 
galaxies. Galaxies are beautiful both in photographs and in VLA 
maps, and they are dynamical objects. I will do this by describing 
an investigation that has been occupying our attention for the past 
couple of years. It is a good example of the method of numerical ex- 
perimentation. 

Experiments in terrestrial laboratories are out of the question 
for problems in the dynamics of galaxies as well. Distances are so 
great that there is no hope of going out beyond the object to get a 
look from the other side. There is no way to kick a galaxy to see if 
it bounces-to check its stability as you would check the stability (or 
robustness) of everyday objects in the laboratory. Even if we could, 
the time scales are so long we could not wait. Dynamical times are 
typically about 100 million years. The need for experimental checks 
on analytical theories is as great in galaxy dynamics as in any area of 
physics. Numerical experiments, based on large scale computations, 
are the closest thing we have to the laboratory experiments of 
other parts of physics. They are much closer in style and spirit 
to laboratory experiments than to conventional (analytic) theory. 
Numerical experiments are initial value problems. They are defined 
by the initial conditions. The same rules are used for carrying out all 
experiments, and those rules contain as much of the essential physics 
of galaxies as possible. Once an experiment is started, it runs to 
“completion” with no interference. 

Studies in the dynamics of galaxies produce strikingly beautiful 
objects as well as charming dynamical insights. Galaxy dynamics is 
challenging because one’s guesses are so often wrong. This happens, 
of course, because self-consistent systems have so many options open 
to them that they usually select one you hadn’t thought about. It is 
exciting as you begin to see which option the system selected and as 


you try to figure out why it chose that option. 

The investigation reported here started out as a study of the 
dynamics of dark lanes in the so-called dark-lane elliptical galaxies. 
This investigation was started by Dr. Althea Wilkinson of the 
University of Manchester in England, and it was carried out by a 
team from which Dr. Bruce F. Smith of NAS A- Ames was a driving 
force. Some totally unexpected and surprising phenomena turned 
up-phenomena that remind one of the strange things observed 
at the center of our own Galaxy. I want to concentrate on those 
phenomena. 

The experiments started with a disk of particles embedded 
within an oblate spheroidal galaxy. Both the disk and the oblate 
spheroid were formed of particles, 400 000 in all. The experiments 
consist of following the dynamical development of this system for 
some time, usually several dynamical time periods, by means of 
a fully self-consistent fully three-dimensional n-body integration. 

The bulk of the mass is in the oblate spheroid, and it represents 
the elliptical galaxy. The disk represents the dark lane. It contains 
little mass (1000 particles), and notwithstanding its dominant optical 
appearance it has little, if any, dynamical effect on the main galaxy. 

The disk was flat at the start. Its normal formed an angle of 45° 
with the spheroid axis. Its center coincided with the center of the 
galaxy. Disk particles orbit around the galaxy center so the disk is 
centrifugally supported against the gravitational forces of the galaxy. 
Disk particles have different velocities at different distances from the 
center, so the disk is said to “rotate differentially.” It also precesses 
differentially, and it soon becomes warped, taking on a sequence of 
beautiful shapes. 

We routinely make motion pictures that show the dynamical 
development of the systems we study. Shapes and motions Eire so 
complex that one cannot understand them otherwise. It is safe to 
say that every significant discovery we have made in 10 years of 
experimentation has started from some feature first noticed in a 
motion picture, whether on film or on a graphics display. Good 
graphics is essential for an experimenter to understand his own 
results. Of course, motion pictures are great for showing results to 
others as well, and a motion picture derived from an experiment run 
on the NAS Cray-2 will be shown. 

The unexpected feature showed up as we were studying an earlier 
version of this system on the then new IRIS graphics displays at 
Ames. We had zoomed in for a tight view of the center, watching 
temporal developments. The center was whipping about. We were 
concerned lest this be a numerical instability. A motion picture 
showing a close-up view of the disk center was made as one of our 
checks. That motion picture will be shown. The motions that 
troubled us axe quite small. The centermost particle moves only 
as far from the center as the radius of a sphere that would contain 
but 4 x 10 -4 of the galaxy’s mass. Nonetheless, it is troubling to 
an experimenter until he can figure out what causes it. If caused 
by a numerical instability, it could signal a trouble with the entire 
sequence of experiments we’ve done over the past 10 years. 

Checks whether we were dealing with a numerical instability 
led us to study / = 1 oscillations of galaxy models, to study the 
stability of the disk with the galaxy immobilized (i.e., the underlying 
galaxy was no longer self-consistent), to study the disk in a harmonic 
oscillator potential rather than the galaxy potential, and to carry 
out a variety of other checks. The / = 1 mode had too low a Q 
to fit the experiment. Experiments with the galaxy immobilized or 
with a harmonic potential showed that self- gravitation of the disk 
was not responsible. Next, we conducted experiments on the galaxy 
alone, with the disk removed. We had to invent a way to look for 
motions of the centred parts in order to do this. The method used 
was to locate the extremum of the gravitational potential. We found 
that the center, as defined by the potential, moved about the center 
as defined by the mass (the centroid). The most dense part of the 
galaxy, indicated by the potential center, need not remain fixed, even 
though the equations of motion require that the mass centroid remain 
fixed. 

Fairly early in this sequence of checks, we plotted the trajectory 
of the normal to a little patch at the center of the disk. That plot 
is shown in Figure 1. It looks like the plots of combined precession 
and nutation in rigid body motion in your elementary mechanics 
text. Indeed, that is just what it shows-combined precession and 
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nutation. It is not as clean as the motion shown in the mechanics 
text, but it is the trajectory of a genuine mechanical motion. It is not 
the trajectory of a numerical instability. 



Figure 1. Trajectory of normal to the center of disk as seen on the 
unit sphere. The trajectory starts at the right side (small amplitude) 
and winds downward and to the left, ending on the left side with 
large amplitude. 


So it’s not a numerical instability. But then, what is it? We 
finally nailed this down with the experiments on the NAS Cray-2. 

The potential near the center of a galaxy is harmonic-it is the 
potential of a (possibly anisotropic) harmonic oscillator. A particle, 
initially at rest at the center, feels tugged first in one direction and 
then in another, by the wanderings of the potential center. It acts 
like a harmonically bound particle in Brownian motion. That is a 
well-known problem. We checked whether this might be the character 
of the motions we were observing by integrating the response of 
a harmonically bound particle to forces it would feel due to the 
wandering center, and found those motions to match the wanderings 
of the centermost particle of our disk fairly well. 

Now we come to the nailing. Experiments with different numbers 
of particles showed that the amplitude by which the potential center 
wanders about the mass center varied inversely with the square 
root of the number of particles in the galaxy. Separate checks, with 
a sequence of initial loads produced by different runs of random 
numbers produced a scatter of potential center positions about the 
same as the early stages of an actual integration. That confirms 
suspicions that y/n noise causes the potential center to differ from 
the mass center. Both checks were run with 25,600, with 100,000, 
and finally with 400,000 particles. The scatter for the integration 
increases at later times. 


The second part of the check is shown in Figure 2. 



Figure 2. Plots of position va. time. Time in the integration is 
plotted horizontally, the x-component of position vertically. The 
motion of the centermost particle is shown in the top trace, the 
potential center wandering in the middle trace, and the motion of a 
harmonically bound test particle subject to forces produced by that 
potential center wandering in the bottom track. All three tracks have 
the same vertical and horizontal scales. 

The actual motion of the centermost particle is shown in the 
topmost trace, and the potential center wandering in the middle 
trace. Finally, a separate integration of a harmonically bound test 
particle under the tuggings of the actual potential center wanderings 
is shown as the bottom trace. Notice the similarity of the top and 
bottom tracks-even to shapes of individual wiggles. The amplitudes 
are a bit different, indicating that something or other about the 
galaxy affects the centermost particle at late times. Further, we 
used oscillation frequencies from an earlier experiment for this 
comparison, rather than re- determining the frequencies from the 
actual experiment. That accounts for the phase difference in the 
two tracks at the end of the experiment. The increasing amplitude 
looks more like a linear growth rate (up to a limiting amplitude) than 
like the square-root of time amplitude of the classical harmonically 
bound particle, suggesting that something more is going on. That 
can be seen in the center track, where one can imagine that the 
potential center wanderings near the end are beginning to follow the 
particle. The system is actually unstable, with feedback coming from 
the potential center’s following the driven particles. But that is a 
physical, not a numerical, instability. 

The experiments are vindicated, but our next task is to deter- 
mine whether this phenomenon is important in real galaxies. As 
mentioned earlier, the wanderings are very suggestive of things that 
happen at the center of our Galaxy, so we could be dealing with an 
important physical efTect. But that involves a separate set of argu- 
ments, one that necessarily hinges on of astronomical details. We 
postpone that to a conference on the dynamics of galactic nuclei. 

I wish to thank the organizers of the NAS Dedication for the 
invitation to participate and to show you some of the beautiful 
results Bruce Smith and I have obtained already with the NAS 
facilities. You can imagine our hopes and excitement over future 
prospects. This work has been partially supported by Cooperative 
Agreement NCC 2-265 between NASA-Ames Research Center and 
the I* diversity of Chicago. 
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